import holoviews as hv
hv.extension('bokeh')
import numpy as np
import scipy.signal
from IPython.display import YouTubeVideo, Audio
1. Introducción al procesamiento digital de señales¶
1.1. Sistemas y señales¶
Definición de sistema: Es un proceso que produce una señal de salida en función de una señal de entrada
Definición de señal: Es una función que describe o representa el comportamiento de un fenómeno físico o sistema. Este comportamiento puede variar ya sea en el tiempo, espacio u otra variable
1.1.1. Ejemplos de señales¶
Revisemos los siguientes ejemplos ¿Cuáles son las variables independientes y dependientes en cada caso?
Las variaciones del IPSA entre 2006 y 2015
Los pulsos de reloj de un microcontrolador
El brillo relativo de una estrella medido desde la Tierra. La estrella es orbitada por un exoplaneta
Una conversación humana grabada digitalmente
La actividad cerebral de una persona medida usando ElectroEncefalografía
Datos meteorológicos de un huracán
Una película
1.1.2. ¿Qué es el procesamiento de señales?¶
Es la disciplina que se dedica al diseño de sistemas para
representar …
filtrar …
codificar …
transmitir …
estimar …
detectar …
inferir …
descubrir …
reconocer …
síntesizar …
reproducir …
etc
… señales
Fuente: Jose S.F. Moura, “What is signal processing”, IEEE signal processing magazine
YouTubeVideo('R90ciUoxcJU')
1.2. Definición matemática de una señal¶
Una señal se representa matematicamente como una función
es decir un mapeo entre el espacio de entrada y el espacio de salida
donde:
\(x\) se llama variable independiente, entrada o argumento y su espacio se llama dominio
\(y\) se llama variable dependiente, salida o retorno y su espacio se llama codominio
1.2.1. Funciones matemáticas típicas¶
Algunas funciones típicas que veremos a lo largo de este curso son
Impulso unitario:
Escalón unitario:
Propiedad: \(\delta[n] = u[n] - u[n-1]\)
Exponencial real:
Exponencial compleja:
Gaussiana
x = np.linspace(-4, 4, num=1000)
p1 = hv.Overlay([hv.Curve((x, x), label='y = x'), hv.Curve((x, np.abs(x)), label='y = |x|')])
p2 = hv.Overlay([hv.Curve((x[x>0], np.sqrt(x[x>0])), label='y = sqrt(x)'),
hv.Curve((x[x>0], np.log(x[x>0])), label='y = log(x)')])
p3 = hv.Overlay([hv.Curve((x, np.exp(x/4)), label='y = exp(x/4)'),
hv.Curve((x, np.exp(-x/4)), label='y = exp(-x/4)')])
p4 = hv.Overlay([hv.Curve((x, np.cos(x)), label='y = cos(x)'),
hv.Curve((x, np.sin(x)), label='y = sin(x)')])
p5 = hv.Overlay([hv.Curve((x, np.heaviside(x, 0)), label='y = step(x)'),
hv.Curve((x, 1.0/(1.0+np.exp(-x))), label='y = logit(x)')])
p6 = hv.Overlay([hv.Curve((x, scipy.signal.unit_impulse(len(x), 'mid')), label='y = delta(x)'),
hv.Curve((x, np.exp(-0.5*x**2)), label='y=exp(-x^2)')])
plot = hv.Layout([p1 + p2 + p3] + [p4 + p5 + p6])
plot.opts(hv.opts.Curve(width=280, height=280, xlim=(-4,4), ylim=(-4,4), show_grid=True, tools=['hover']),
hv.opts.Overlay(legend_position='top')).cols(3)
1.2.2. Algunas propiedades de las funciones¶
Función par o simétrica:
Función inpar o antisimétrica:
Función periódica: Sea \(T\) el período de la función
Función lineal: Sean \(a\), \(b\) coeficientes escalares
El periódo de una función es equivalente al inverso de su frecuencia fundamental
La siguiente animación muestra como varía una sinusoide en función de su frecuencia
time = np.linspace(-6, 6, num=1000)
hmap = hv.HoloMap(kdims=['Frecuencia'])
for freq in np.arange(0, 1, step=0.01):
hmap[freq] = hv.Curve((time, np.sin(2*np.pi*freq*time)))
hmap.opts(hv.opts.Curve(width=500, height=300))
hv.output(hmap, holomap='gif', fps=10)
1.3. Señales analógicas y digitales¶
La siguiente figura muestra cuatro combinaciones en que se pueden presentar las variables independiente y dependiente de una señal
Variable independiente continua (columna izquierda)
Variable independiente discreta (columna derecha)
Variable dependiente continua o análoga (fila superior)
Variable dependiente discreta o cuantizada (fila inferior)
Las señales “naturales” son en general analógicas de tiempo continuo. Sin embargo los computadores trabajan con señales digitales de tiempo discreto
Podemos convertir una señal analógica a una representación digital con los siguientes pasos:
Se discretiza en el tiempo muestreando según el tiempo de muestreo \(T_s\) del sistema, por ejemplo:
Esto convierte la señal analógica \(x_A\) en una secuencia discreta de valores \(x\) indexada por un entero \(k\), es decir un arreglo
Se cuantiza la amplitud de la señal, por ejemplo:
Esto convierte la variable dependiente continua (en este caso voltaje) a una representación discreta (en este caso binaria)
Ts = 0.249
x = np.linspace(0.0, 2.0, num=1000)
y = np.sin(2.0*np.pi*1.0*x)
p1 = hv.Curve((x, y), label='Señal analógica').opts(width=250)
p2 = hv.Curve((x, Ts*np.round(y//Ts)), label='Señal cuantizada').opts(width=250)
p3 = hv.Points((x[::20], Ts*np.round(y[::20]//Ts)), label='Señal digital').opts(width=250)
(p1 + p2 + p3)
1.4. Propiedades fundamentales de las señales¶
En general trabajaremos con señales que, dentro de un rango de interés, son acotadas en energía, potencia y/o en ancho de banda. Es decir son señales finitas que no divergen ni se vuelven singulares
La energía de una señal mide su “tamaño” o el “espacio que ocupa”. Para una señal analógica y discreta, respectivamente
Matemáticamente, se define como
Una señal acotada en energía debe cumplir \(E_s < \infty\)
La potencia promedio de una señal se define como su energía por unidad de tiempo
Matemáticamente, se define como
Una señal acotada en potencia debe cumplir \(P_s < \infty\) y una señal de energía finita tiene potencia cero
Por otro lado una señal de potencia finita tiene duración infinita
El ancho de banda de una señal mide su tasa de cambio o velocidad
Una señal acotada en ancho de banda debe tener transiciones suaves
Ejemplo formativo
¿Cual es la energía de esta señal?
dt = 0.001
x = np.arange(-5, 20, step=dt)
y = np.zeros(shape=(x.shape))
y[x>=0] = 2*np.exp(-x[x>=0]/2)
# Integral numérica
print((y**2).sum()*dt)
hv.Curve((x, y)).opts(width=500)
4.0020003250765726
¿Puedes comprobar el resultado resolviendo la integral?
1.5. Procesamiento estadístico de señales¶
Las funciones matemáticas que vimos anteriormente son ejemplos de señales deterministas. Una señal es determinista si puede describirse completamente por una ecuación matemática
Por ejemplo
x = np.linspace(-10, 10, num=1000)
y = np.cos(x)*np.sin(2*x) + 0.1*x
hv.Curve((x, y)).opts(width=500)
Sin embargo, existen algunos fenómenos que son difíciles o imposibles de describir de forma determinísta, como son por ejemplo las señales estocásticas
Una señal estocástica tiene un cierto grado de aleatoridad, para describirla debemos usar teoría de probabilidades
Por ejemplo
quiere decir que \(z\) se distribuye normal con media \(\mu\) y covarianza \(\Sigma\)
Podemos estudiar la distribución de \(z\) usando un histograma
z = np.random.multivariate_normal(np.zeros(100), np.eye(100))
edges, bins = np.histogram(z, bins=20)
hv.Histogram((edges, bins)).opts(width=500)
1.5.1. Breve repaso de estadística¶
Variable aleatoria (VA): Es una variable cuyos valores posibles son resultados de un fenomeno aleatorio
Se describe en términos de su distribución
función de densidad de probabilidad (fdp) para variables continuas
función de masa de probabilidad (fmp) para variables discretas
Los valores observados a partir de la VA se llaman realizaciones
Por lo general asumimos que las realizaciones son iid: independientes e identicamente distribuidas
Usualmente se denota como \(X\) (mayúscula) y sus realizaciones como \(\{x_1, x_2, \ldots, x_N\}\) (minúscula)
Una distribución se describe a través de sus momentos estadísticos. Para una variable aleatoria \(X\) su momento de orden \(k\) es
y su momento central de orden \(k\)
donde la esperanza se define como
y la media muestreal es
Los momentos estadísticos más usados son la
Media (\(\mu_1\)): Define la localización de la distribución
Varianza (\(\hat \mu_2\)): Define la dispersión o ancho de la distribución
Simetría (\(\frac{\hat \mu_3 }{\sqrt{\hat \mu_2^3}}\)): Define el peso relativo de las colas de la distribución
Kurtosis (\(\frac{\hat \mu_4 }{ \hat \mu_2^2}\)): Define el peso relativo entre las colas y el centro (moda) de la distribución
La siguiente gráfica muestra como cambia una distribución ante estos momentos
1.5.2. Proceso aleatorio/estocástico¶
Un proceso aleatorio (PA) es una colección de variables aleatorias indexadas
Llamamos a una realización del PA: serie de tiempo
Un PA se dice estacionario si sus momentos no varían con el tiempo, es decir
Un PA se dice ergódico si los momentos se pueden deducir a partir de una realización (suficientemente larga) del mismo
Esto es útil pues muchas veces observamos sólo una realización
Los promedios muestreales del ensamble equivalen a promedios muestreales temporales
La siguiente figura muestra la diferencia entre un promedio temporal y un promedio de las realizaciones
Ejemplo: Realizaciones de un PA con distribución Gaussiana
Mientras más muestras observamos, mejor se aprecia la distribución subyacente
np.random.seed(12345)
data = np.random.randn(10000)
x = np.linspace(-3, 3, num=100);
signal = hv.HoloMap(kdims=['Muestras'])
histogram = hv.HoloMap(kdims=['Muestras'])
for N in [10, 50, 100, 500, 1000, 5000, 10000]:
signal[N] = hv.Curve((data[:N]), kdims='tiempo', vdims='datos').opts(xlim=(1, N+1), logx=True)
edges, bins = np.histogram(data[:N], bins=20, range=(-3, 3), density=True)
histogram[N] = hv.Overlay([hv.Histogram((edges, bins), kdims='datos', vdims='Frecuencia').opts(alpha=0.5),
hv.Curve((x, np.exp(-x**2)/np.sqrt(2.0*np.pi)))])
hv.output((signal + histogram), holomap='gif', fps=1)
Ejemplo: Realizaciones de un PA con distribución Normal/Gaussiana con media cero
Cada columna tiene una covarianza distinta. Mientras más grande es el valor de “tau” más correlacionadass están las muestras
np.random.seed(12345)
dt = np.repeat(np.reshape(np.arange(100), (1, -1)), 100, axis=0)
def draw_random_plot(tau):
p = []
for k in range(0, 25, 5):
Sigma = np.exp(-0.5*np.square(dt - dt.T)/tau**2)
data = k + np.random.multivariate_normal(np.zeros(100), Sigma)
p.append(hv.Curve((data)))
return hv.Overlay(p, label=f'tau = {tau}').opts(width=250, height=400, yaxis=None, xaxis=None)
draw_random_plot(0.5) + draw_random_plot(1) + draw_random_plot(2)
1.5.3. Ruido¶
En la práctica no solemos observar señales puramente deterministas. Una de las causas de esto es la presencia del ruido
El ruido es una corrupción indeseable y usualmente estocástica que modifica la señal de interés
Corrupción aditiva
Un supuesto muy utilizado es considerar que el ruido es una corrupción aditiva
donde \(y\) es la señal observada, \(x\) es la señal de interés y \(n\) es una señal de ruido
np.random.seed(12345)
t = np.linspace(-10, 10, num=500)
x = np.cos(t)*np.sin(2*t) + 0.1*t
s = 0.5
n = s*np.random.randn(len(t))
print(f"SNR: {10*np.log10(np.sum(x**2)/np.sum(n**2))}")
p1 = hv.Curve((t, x))
p2 = hv.Curve((t, n))
p3 = hv.Curve((t, x+n))
(p1+p2+p3).opts(hv.opts.Curve(width=250, height=250))
SNR: 4.063138926337416
Podemos cuantificar el “nivel de ruido” o la claridad de la señal observado en términos de la razón señal a ruido (Signal to Noise Ratio, SNR)
La SNR se define como la razón entre la energía de la señal y la energía del ruido
La SNR se mide en decibeles [dB]
1.5.4. Tipos de ruido¶
Podemos clasificar el ruido según
la distribución que sigue, por ejemplo ruido Gaussiano o ruido Uniforme
la dependencia temporal/espectral entre sus realizaciones, por ejemplo ruido blanco o ruido rojo
¿Cómo se ve y suena un ruido blanco?
Fs = 44100
y = np.random.randn(Fs*5)
y = y/np.max(np.absolute(y))
hv.Curve((y), 'tiempo', 'ruido').opts(width=500)
¿y un ruido rojo?
dy = np.random.uniform(-1, 1, Fs*5)
y = np.cumsum(dy)
y = (y - np.mean(y))
y = y/np.amax(np.absolute(y))
hv.Curve((y), 'tiempo', 'ruido').opts(width=500)
Audio(y, rate=Fs)